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Abstract 

A weak wave turbulence theory is developed for two-dimensional (2D) magnetohydrodynamics 
(MHD). We derive and analyze the kinetic equation describing the three- wave interactions of 
pseudo-Alfven waves. Our analysis is greatly helped by the fortunate fact that in 2D the wave- 
kinetic equation is integrable. In contrast with the 3D case, in 2D the wave interactions are 
nonlocal. Another distinct feature is that strong derivatives of spectra tend to appear in the 
region of small parallel (i.e. along the uniform magnetic field direction) wavenumbers leading to 
a breakdown of the weak turbulence description in this region. We develop a qualitative theory 
beyond weak turbulence describing subsequent evolution and formation of a steady state. 
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I. INTRODUCTION 



Magneto-Hydro-Dynamics (MHD) is of great interest for modeling turbulence in magneti- 
cally confined and unconfined plasmas. In astrophysics its applications range from solar wind 

to the Sun {2], to the interstellar medium and beyond Q. At the same time, MHD is 
also relevant to large-scale motion in nuclear fusion devices such as tokamaks lal. One of the 

n 

pioneering results of incompressible MHD turbulence has been obtained by Iroshnikov |6[ and 
Kraichnan [3] (thereafter IK) who proposed an extension of the Kolmogorov phenomenology 
8], originally derived for Navier-Stokes equations for Hydro-Dynamics (HD). For simplicity 
the assumptions of homogeneity and isotropy were made by Kolmogorov, and the energy 
cascade was supposed to be dominated by local (in scale) interactions between eddies of 
similar size. Then the Kolmogorov phenomenology leads to the well-known one-dimensional 
kinetic energy spectrum E(k) ~ k~ 5 ^ 3 , where k is the wave vector. The associated cascade 
properties for its inviscid invariants differ for 3D and 2D turbulence: while in 3D the energy 
and the kinetic helicity exhibit direct cascades, in 2D the energy cascades inversely - still 
with a —5/3 scaling - whereas a direct cascade is found for the enstrophy which leads to 
spectrum E{k) ~ k~ 3 at the small scales. 

IK modified the Kolmogorov phenomenology by taking into account the magnetic field. 
They also assumed turbulence homogeneity, isotropy and locality of interactions. However, 
there exist fundamental differences between the Kolmogorov and the IK theories. First of 
all, in MHD the energy cascade is supposed to be dominated not by the interactions between 
eddies but between Alfven wave packets propagating in opposite directions; this modification 
leads to the energy spectrum E{k) ~ k~ 3 / 2 . Furthermore, unlike hydrodynamics the cascades 
of the ideal MHD invariants exhibit some similarities in their behaviour in 2D and 3D 
in both cases, the cascades have the same direction with a direct cascade for the energy and 
cross-helicity, and an inverse cascade for the magnetic helicity or the anastrophy. 

The differences between MHD and HD turbulence go beyond these classical properties. 
In the IK theory the large-scale magnetic field is supposed to play the role of external field 
which is necessary for the existence of Alfven waves but its main effect, i.e. anisotropy, is not 
taken into account. The importance of an external magnetic field has been discussed many 



times during the last two decades 1 OMl 81] and the anisotropic behavior has been shown in 



direct numerical simulations for both 2D 



Uj and 3D [12|. 



Despite some similarities - like the cascade directions - the question about the identifica- 
tion of differences between 2D and 3D MHD turbulence still represents an important issue. 



In early numerical studies 11] mainly 2D simulations were performed because of the limited 



numerical resources available and the non-accessibility of 3D calculations at high Reynolds 



numbers. Nowadays, 3D MHD numerical simulations are commonly achieved 26 1 . At 
;he same time, 2D simulations are still used for the illustration of new numerical techniques 
27j. There is also an interest for the investigation of freely decaying MHD turbulence be- 



cause the Reynolds numbers can be higher in 2D than in 3D [23|. When the magnetic 
field perturbations are small compare to a uniform background magnetic field, the 2D MHD 
equations are sometimes used to model turbulence. Such a situation is particularly relevant 
for solar coronal loops as well as for reduced models for fusion plasmas in slab geometry jlj] . 

Strong statements were made by some authors that 2D simulations can be safely used 
to model 3D situations because the properties of the 2D and the 3D MHD turbulence are 



essentially the same 23_|, |28| . One of the motivations of the present paper is to test validity 
of this claim in a special case when the external magnetic field is strong. 

In this paper, we consider 2D MHD in the presence of a strong background magnetic field 
which implies realization of the weak turbulence regime. One of the main advantages of this 
regime is the fact that it allows one to derive accurate analytical results for the spectrum. An 
explicit comparison will be made between the weak turbulence regimes in 2D and in 3D; the 



latter was analyzed rigorously in [29 L The weak wave turbulence approach is widely familiar 
to the plasma physics community 3Ch42||. It is a statistical description of a large ensemble 
of weakly interacting dispersive waves. The formalism leads to wave kinetic equations from 
which exact power law solutions can be found for the energy spectra. There were several 
reasons which postponed development of weak turbulence for Alfven waves. The first one is 
their semi- dispersive nature. Typically, the wave-kinetic approach cannot be used for non- 
dispersive waves since such wavepackets propagate with the same group velocity even if their 
wavenumbers are different, the energy exchange between such waves may not be considered 
small, and may lead to possible energy accumulation over a long time of interaction. The 
Alfven waves represent a unique exception to this rule because co-propagating waves do 
not interact and the nonlinear interaction is present only for counter-propagating waves. 
The latter pass through each other in some finite time and no long time cumulative effect 
occur. That is why the Alfven waves represent a unique example of semi- dispersive waves 



for which the wave turbulence theory applies. The second reason which renders the weak 
turbulence theory for Alfven waves very subtle is the fact that the domination of three- 
wave interactions - as assumed by IK - may be questionable. While in [43J the three- wave 



interactions were declared absent, the IK argument has been re-established i n |13l . |29 



The weak turbulence theory for 3D incompressible MHD was developed in [29[ (see also 



45 



46| ) where three- wave kinetic equations were derived with their exact solutions via a 



systematic asymptotic expansion in powers of small nonlinearities. 

The main goal of the present paper is to derive the weak turbulence equation for the 
2D MHD, analyze it and make a comparison with the 3D case. The crucial technical step 
which allows a comprehensive theoretical analysis of the solutions consists in transforming 
the wave kinetic equation into an integrable form by Fourier transforming it and separating 
the transverse and the parallel dynamics by using a self-similar "effective time" variable. 

This article is organized as follows. In section [Til we derive the weak turbulence kinetic 
equation through a general perturbative procedure. In sections [TTTJ and HV] we proceed with a 
detailed investigation of its properties in an anisotropic limit. The main goal of such a study 
is to verify whether or not the turbulence is local. First of all, we will consider the steady 
state behavior by looking for Kolmogorov-Zakharov type solutions and check their locality. 
Next, we will proceed with investigation of an unsteady spectrum evolution by considering 
two different cases with a gaussian-shaped source and different kinds of dissipation: a uniform 
friction and a viscosity. Due to integrability of the weak wave-kinetic equation in the first 
case it is possible to find an exact solution. In the second case, a qualitative analysis for 
the steady state is complemented by a numerical simulation of the spectrum evolution. The 
goal of section [V] is to develop some qualitative reasoning about the turbulent behavior of 
our system near the applicability margin of wave-kinetic formalism and beyond. Formation 
of steady state is also discussed. Finally, we present a summary of our results in section IVII 



II. WAVE-KINETIC DESCRIPTION 



A. Alfven waves 



In 3D incompressible MHD there exist two different kinds of Alfven waves 15j. The 



first kind, called Shear-Alfven waves (SAW), have fluctuations of velocity and magnetic field 



transverse to the initial background magnetic field Bo , whilst the other kind, called Pseudo- 
Alfven waves (PAW), have fluctuations along B . Both waves propagate along B at the 
same Alfven speed. 

The weak wave turbulence approach for incompressible MHD applies for a small nonlin- 
earity, e ~ b±k±/B k\\ <C 1, were b± is the perpendicular magnetic field perturbation and 
fcn and k± are the wave vector components in the parallel and the perpendicular directions 
to Bo- Additionally, a strong anisotropy condition is often used, a = k±/ku ^> 1. In the 
3D case, it was shown that in the leading order of weak nonlinearity, e < 1, and strong 
anisotropy, a 3> 1, the SAW interact only among themselves and evolve independently 
from the PAW. At the same time, the PAW scatter from the SAW without amplification or 
damping, and they do not interact with each other. Such a behaviour does not rule out a 
possibility for the PAW to be interacting among themselves in the next order of expansion 
in \ jo. However in the 3D case, such a process is sub-dominant to a stronger interaction 
with the SAW and was not considered yet. 

In the 2D case, due to the geometrical restrictions, it is only possible to have the PAW 
and not SAW. In this paper we will see that three-wave interactions of PAW do occur in the 
2D case in the next order of expansion in 1 jo and represent the dominant process in the 
nonlinear evolution. 

B. Interaction representation. 

The ideal incompressible MHD system in Elsasser variables z = v + sb, with s = ±1, is 
given by [28( 



where v is the fluid velocity, b is the magnetic field fluctuation (in velocity units), B is a 
uniform background magnetic field (also in velocity units, i.e. the Alfven speed) and P* is 
the total (thermal plus magnetic) pressure. In what follows we suppose that the background 
magnetic field is directed along the x axis, B = B x. In the coordinate notations we have 



(d t - sB ■ V + z- s ■ V) z s = -VP 



(1) 
(2) 



V -z s 







(d t - sB d x ) z) 



£ Z n O n Zj OjP*, 



(3) 



The nonlinear terms in eq.([3]) include Elsasser variables of the opposite signs only. Therefore, 
the nonlinear interactions take place only between counter-propagating waves. 

The first step in the general procedure of the wave kinetic formalism is to identify the 
linear modes. Neglecting the nonlinear terms in the right hand side (r.h.s.) of the equation 
(J2J) (which includes the pressure term) and looking for solutions in the form of a wave, 

Z s . ~ e i{k x x+k y y)-iu s t ^ 

we get two linear modes, 

uj s = —sB k x , s = ±1, (6) 

which propagate parallel to the background magnetic field (in both directions) with group 
velocity = — sB . Let us suppose that our system is periodic in the physical space with 
period L in both x and y, and let us introduce the Fourier series: 



(x,t) = ^ S s (k,t)e 4k - x , (7) 

k 

where wave vector k takes values on a 2D grid, k = (k x ,k y ) = (2irm x / L,2irm y / L) where 
m x , m y G Z. Then, by applying the divergence operation on both sides of the equation 
and by using fll]), we find the expression for the Fourier coefficients of pressure P*: 

A(k) = -AT 2 ( k 2 • a ~ S ( k i> *))( k ' aS ( k 2> *))<Kki + k 2 - k), (8) 

ki,k 2 

where S(p) is the Kronecker delta: 5(p) = 1 for p = and zero otherwise. Thus, equation 
03]) in Fourier space becomes: 



id t - u s ) a s (k, t) = ( k • a ~ s ( k i> *)) 

ki,k2 



a s (k 2 ,t)-^(k-a s (k 2 ,t)) 



S(ki + k 2 -k). (9) 



Using the incompressibility condition 



we reduce (jSJ) to one scalar equation, 

(id t -u s )al(k,t) = eJ2 J k y (kX fei kl ^ (k fe2 k2) ^(ka^jdkx dk 2 . (11) 



Let us now introduce the representation of interaction variables, 

. . k 



c k = i 



ky€ 



a x (k,t)e^\ (12) 



which represent slowly varying wave amplitudes. Factor e lujSt here compensates for the fast- 
scale oscillations arising due to the linear dynamics. We have introduced a formal constant 
small parameter e C 1 for easier counting of the powers of nonlinearity, assuming now that 
Ct ~ 1. Then the system in the interaction representation variables is: 

dtct = e V ™ c t 4/ iklxt ^i + k 2 - k), (13) 

ki,k2 

with the interaction coefficient: 

_ (k-k 2 )[k 1 xk 2 ] z 
Vl2k ~ khk 2 • (14) 

Note that so far we have not used smallness of e and the eq. (I13p is completely equivalent to 
the initial system ©-(jll)- 

C. Wave- kinetic equation 



The standard weak turbulence approach 35|, |46( exploits the smallness of nonlinearity, 



randomness of phases and the infinite box limit. In Appendix [AJ we are applying this 
approach to the system (I13p and obtain the following kinetic equation for the wave spectrum 



* J Vfe 2 12 < [ n ± - n£] 8 (k - kx - k 2 ) 5(2k lx )dk x dk 2 , (15) 

where the interaction coefficient is as in expression ( fl4l) . Here, S(p) means Dirac's delta 
function. In the following sections we proceed with detailed analysis of this equation. 

D. Anisotropic limit 

One remarkable property of MHD turbulence, which makes it very different from the 
HD one, is its strong anisotropy in presence of strong background mag netic field. It was 
illustrated with direct numerical simulations in both 2D 11] and 3D [12j. The wave-kinetic 



formalism confirms such an anisotropy through the kinetic equation structure. In fact for 



Alfven waves, the resonant three-wave interaction [11[ is organised in such a way that one 
member of each triad must have its wave vector perpendicular to the external magnetic 
field. At the same time, the two other waves in the same triad must have their parallel 
wavenumbers equal to each other: fen = k 2 \\. Formally, this is seen in both the 2D and 



the 3D kinetic equations whose r.h.s. contains a delta function 8(2ki\\), see equation ([15 



for the 2D case and the equation (26) in 29( for the 3D case. Using this delta function 



and integrating over feiii we see that the parallel component of the wavenumber enters into 
the kinetic equation as an external parameter and the spectrum dynamics is decoupled at 
each level of ku. In other words, there is no energy transfer in the parallel (to the external 
field Bo) direction in the k-space. The initial spectrum is spreading over the transverse 
wavenumbers k±, and not in the k\\ direction. For large time, such a spectrum becomes very 
flat, pancake-like. 

The two-dimensionalisation of the total energy means that, for large time, the energy 
spectrum is supported on a volume of wave- numbers such that for most of them k± 3> fc||. 

Thus, let us consider the anisotropic limit for kinetic equation ( I15p . which reads in 2D 
as k y 2> k x . Taking into the account the resonant interaction conditions for the parallel 
wavenumbers, we will have a considerable simplification of the interaction coefficient: 

V k i2 = K, (16) 

and the kinetic equation will take the following form, 

^(kx, k y ) = Tik 2 x I n T (0, h y )[n ± (k x , k 2y ) - n^ftx, k y )}5(k y - k ly - k 2 y)dkiydk 2 y. (17) 



This equation describes three- wave interactions of PAW in 2D in the anisotropic limit. One 
can immediately see that the energy is conserved separately in the "+" and "-" waves 
separately at each k x : 

d t J n T (k x ,k y ,t)dky = 0. (18) 

Factor k 2 x in the r.h.s. of relation (ITT)) is very important. In the 3D case, there exists 
a similar term which corresponds to a sub-leading contribution. We remind that in 3D in 
the leading order of the perturbation theory, the PAW are scattered on SAW and do not 
interact directly to each other. In the 2D, there is no SAW and, therefore, the r.h.s. of ( II 7p 
becomes the leading order contribution. 

Further, in leading order of 3D there is no k\ factor, and substitution n(k±,k\\,t) = 
n±(k±,t)n\\(k\\) leads to an equation for n±(k±,t) which does not involve k\\. In 2D, one 
can also obtain an equation which does not involve k x but for this, one has to introduce an 

8 



"effective time" variable r = nk x t: 

d T n ± (k xl k y ; r) = J n T (0, k ly ; 0) [rr^A^, k 2y ; r) - n ± (A; ;r , r)] 

x5 (k y - k ly - k 2y ) dki y dk 2y . (19) 

Now we can seek solution in the following form, 

n(k x , k y , r) = fx(k x )r](k y , r), (20) 

where fi(k x ) represents the parallel (non-evolving) component of the energy spectrum and 
i](k y ,T) is the perpendicular one. Without loss of generality we can assume /i(0) = 1. 
Substituting (1211]) into f lT^j) we have the following equation for i], 

d T r] ± (k y , t) = J i] T (k ly , 0) [i] ± (k 2y , r) - ^(ky, r)] S (k y - k ly - k 2y ) dk ly dk 2y . (21) 

Thus, like in 3D, we have an evolution equation for the perpendicular part of the spectrum 
which does not explicitly depend on the parallel part. However, the qualitative difference 
with the 2D case is that there is an implicit dependence on k x via the effective time variable 
r, which in particular leads to the fact that in the r.h.s. of equation f[2"Tj) one of r/'s is taken 
at t = 0, making this equation linear and, as we will soon see, integrable. Another distinct 
feature arising from such an implicit dependence on k x via r is sharpening of the spectrum 
at small k x leading to the breakdown of the wave-kinetic description. Later in this paper we 
will study this effect and its consequences. 



III. KOLMOGOROV-ZAKHAROV SPECTRA AND THEIR LOCALITY 

At the first step of our investigation of the wave- kinetic equation (12 ip. we seek exact 
stationary power law solutions, r)(k y ; oo) ± oc ky . For this, we will use so-called Zakharov 
transformation, 

h h, k 2 

We split the integral into the r.h.s. of (12 ip in two equal parts and we perform the Zakharov 
transformation in the integrand of one of these parts. This leads to the following conditions 
on the exponents, 



-2. (23) 



However, the resulting power law spectra represent true mathematical solutions if and only 
if the original (before the transformation) integral in the r.h.s. of (|2"T1) converges at these 
solutions. In Appendix |Bj we perform such a convergence study and show that this integral 
never converges and therefore, no Kolmogorov-Zakharov solutions are possible for 2D MHD 
system, contrary to the 3D case. We remark that the balanced turbulence spectrum with 
u + = v~ = —1 has a logarithmic divergence. Thus, one could anticipate that such a marginal 
nonlocality could be "fixed" by logarithmic corrections. We will see later that this is false. 



IV. INTEGRATION OF THE KINETIC EQUATION 

A remarkable property of the kinetic equation (1191) is its simplicity. In this section we 
show that in some physical situations it can be solved analytically. Let us introduce into 
the equation ffl9l) sources and sinks of waves: 

d T n ± (k x , k y , r) = J n T (0, k ly , 0) [^(k^ k 2y , r) - n ± (A; a ., k y , r)] 

x8(k y - k X y - k 2y ) dk ly dk 2y + T{k y , k x ) - a d n(k x , k y , r). (24) 

The function F(k x , k y ) may represent forcing or dissipation, depending on the choice of the 
sign before it, and the constant ad introduces a uniform friction. 

In order to use factorisation fTSOj) and eliminate fi(k x ) in both sides of the forced ki- 
netic equation, we assume the following type of force/dissipation function, Tiky^kx) = 
J~x{k x )J~y{ky). Then the parallel component of ( 1201) must be chosen as fi(k x ) = Fxikx). 
Finally, we obtain the following forced/dissipated kinetic equation for the perpendicular 
component of the energy spectrum, 

<9 r ?7 ± (fc 2/ ,r) = J r] T (k ly ,0) [rj ± {k 2y ,T) -^(fcy.r)] 

xS(k y - ki y - k 2y ) dkiydk 2 y + Fy{k y ) - ^(ky, t) . (25) 



A. Pseudo-physical space 

A considerable simplification of equation (J25l) can be obtained with performing the inverse 
Fourier transform on r}{k y: r): 

£ ± (y,r) = J r ] ± (ky,r)e lk ^dk y . (26) 
10 



We will call E ± (y, r) the pseudo-physical space energy, keeping in mind that what is trans- 
formed is the spectrum and not the original wave variable. 

We arrive at the following representation of equation (12"5|) in the pseudo-physical space, 

d T £± (y, r) = £T (y, T ) [£± (y, 0) - E ± (0, 0) - a d \ + F{y). (27) 
B. General solutions 

Let us consider the balanced turbulence case with £ + (y, r) = £~(y, r). Then, the general 
solution of equation (1271) can be written as: 



£(y, r) = C (y)e^-^-^ - ?M , ( 28 ) 

£{y,0) - £ (0,0) - a d 

where the first term represents the general solution for the homogeneous equation and the 
second term is a particular (time independent) solution of the inhomogeneous equation. 
Function C(y) has to be fixed by the initial condition, 

C(y) = £(y, 0) + H ^ . (29) 

£{y,0) - £(0,0) - a d 

Now let us consider two particular examples of the forcing and dissipation. In both 
cases we will assume a Gaussian shape forcing, J-(y) = &fe kfV where constants <jf and kf 
represent the forcing strength and its characteristic wave vector respectively (in the ky-sp&ce 
the forcing is also Gaussian, centered at k y = and with width kf). In the first example, the 
dissipation will be represented by a uniform friction. Here, we can write analytical solutions 
of the kinetic equation on the pseudo-physical space. In the second case, we will consider 
a viscous dissipation. Here, a qualitative analysis of the stationary regime can be done. In 
order to illustrate the spectrum evolution in that case, a numerical solution will be used. 

1. Uniform friction. 

For the uniform friction case we have: 

£ (y, T ) = C (y)e^- £ ^-^ T - ^T, '\ • (30) 

1 ^ ^ ' £(y,0)-£(0,0)-a d 1 ; 

For simplicity, let us use a single-wave initial condition 

£(y,0) = 2Acos (k y) , A = const > 0, (31) 
11 



which corresponds to two (^-functions in fc^-space: at k y = ±&o- 

Then, we can find a function C(y) using equation ( )29|) and substitute it into our solution, 
which yields: 

-kjy 2 /2 



2A cos(k y) + <Tj ' 



2A (cos(k y) - 1) - a d 



-k 2 f y 2 /2 

<jfe i y ' 



(32) 



2 A (cos(k y) - 1) - a d 

Let us examine the steady state, which corresponds to the limit t — > oo and, therefore, 
r — > oo. Note, however, that the time for the steady state to form becomes longer as k x gets 
less, and there are always very small k x where the spectrum is evolving at any large time. 
In the limit r — > oo the solution is given by the second term in the r.h.s. of f )32|) . Far from 
the initial and the forcing scales, at k 3> k$ and k ^> kf, which corresponds to y ^ 1/ko 
and y < l/k f , we have cos(k y) = 1 - (k y) 2 /2 + O{(k y) A ) and e" fe ^ 2/2 = 1 + 0{{k f y) 2 ). 
Thus for this range of scales we have the following expression for the steady state solution 
in the pseudo-Fourier space, 

g ( y oo) = — ^^7, where A = Akl. (33) 

Performing Fourier transform of S(y, oo) we get the steady spectrum, 

1 f°° 

r) (k y , oo) = — / E(y, oo) e- ik ^dy. (34) 



For wavenumbers in the inertial range k ,kf k kd = y/X/od, expression ( 1331) becomes 
effectively a delta function in the integrand of fl3] 



°d + Xy 

and we have 

T] (k v ; oo) = --^L= = - y a/ . (36) 

Therefore we can reach the conclusion that in the equilibrium state the energy spectrum of 
our system in the inertial range is flat. Formally, it is a power law with exponent v = 
which is very different from the Kolmogorov-Zakharov exponent v — — 1 found in section 
IIHI Recall that the Kolmogorov-Zakharov in the balanced case was found to be marginally 
nonlocal and the common wisdom would suggest that it could be fixed by a log correction. 

12 



As we now see this is false: our exact solution has a completely different exponent and has 
no log factor. 

We also see that our exact solution is nonlocal: it does not just depend on the energy 
flux but it contains information about both the sources and the sinks as well as about the 
initial conditions. 



2. Viscous friction. 

Let us now replace the uniform friction by a viscous dissipation keeping the same one- wave 
initial condition as before. Equation ( 127]) becomes: 



dS(y; 



j_ 



_ 2A (cos(fc y) - 1) £(y; r) + o v ^ 1 + a f e~—, (37) 

where a v denotes now the viscosity coefficient and we have used the initial conditions ( l3Tj) . 

To realise this estimation, we need first to get the expression for the steady state solution. 
Let us examine the steady state concentrating, like before in the uniform friction case, on 
the scales less than the forcing and the initial scales, which in terms of the pseudo-physical 
space variables means y -C l/k and y <C l/kf. Performing the same type of expansion in 
small y as before, we have 

d 2 £(y; oo) 



dy 2 

By performing the following rescaling 



\y 2 £{y;oc,)+a, = 0. (38) 



y = y{ ±Y,£ = ^£ (39) 



we obtain 

d 2 £(y, oo) 



y 2 £(y,oc) + l = 0. (40) 

The homogeneous part here is the equation of the parabolic cylinder. Its solutions are the 
parabolic cylinder special functions, whose properties and asymptotics can be found e.g. 
in (47)]. Qualitatively, the behaviour is similar to the one we had in the previous (friction 
dissipation) example: it has a maximum at y = and it decays for y — > 00 (faster than in 
the previous example). In fig. we present such a solution in the pseudo-physical space 
obtained using Matlab in the interval with y 6 [—6.1,6.1]. In order to obtain the decaying 
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y 



FIG. 1. Stationary solution in the pseudo physical space. 



solution we need to take S(0) = 1.3110288959 with high accuracy (to eliminate contribution 
of the growing parabolic cylinder function). 

Respectively, in the fc^-space we again have a flat spectrum in the inertial range, 



where C is an order one constant. Once again we see that the spectrum is nonlocal (i.e. 
it is dependent on the details of the forcing and the sink parameters rather than just the 
energy flux), and its exponent is zero (i.e. it is not a log corrected Kolmogorov-Zakharov 
spectrum). 

In order to illustrate the dynamical evolution of the spectrum in the viscous dissipation 
case, we perform direct numerical simulations for the kinetic equation ( 1251) in the ky-spa.ce. 
In this equation we take i] + = r]~, <r d = and J- y {k y ) = T force — a uky with a u = 10~ 6 . The 
result is presented in Fig. ([2]). The black curve on this figure represents the initial spectrum 
with a Gaussian shape with a large-scale forcing J- 'force = 3 * 10 _4 /fc realised for k e [3,9]. 
We see that the system converges toward a steady state with a flat spectrum in the inertial 
range. 

V. BEYOND WEAK TURBULENCE. 

In this section we are providing, at a qualitative level of rigor, a description for the energy 
spectrum behaviour beyond the weak turbulence regime at the late stage of the evolution. 
For the wave-kinetic equation to be valid the nonlinearity has to remain weak, i.e. the 
nonlinear time scale should be much longer than the linear wave period. This results in the 




for 




(41) 
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FIG. 2. Spectrum dynamics for a viscous dissipation case r = 10000. 
following condition, 



< 1. (42) 



tL B k x 

Here we have used the estimate for the nonlinear evolution time taking the standard hydro- 
dynamic non-linear time scale. In our dynamical equations, the nonlinearity related to the 
magnetic field is of the same order. Finally, the applicability condition can be rewritten as 
a limitation for the parallel wave number: 

k x >K = b -^. (43) 

-DO 

It means that for small values of parallel wave vector component near k* the kinetic equation 
becomes invalid. The question about applicability of the wave- kinetic equation near k\\ = 



has been frequently discussed in the literature, eg. in the 3D MHD turbulence context [29 ]. 



In particular it was speculated in 29] that a sufficient spectrum smoothness for wavenumbers 
near ku = must be present for the wave- kinetic equation to be applicable. In the 2D case, 
the smoothness of spectrum near k x = is asymptotically (in time) broken due to a special 
structure of the kinetic equation. 

Indeed, as we mentioned before, the wave-kinetic equation for the 2D MHD system is 
formulated in terms of a self-similar "time" variable r = k 2 x t. Therefore dependence on the 
parallel wavenumber is still present in the perpendicular part of the energy spectrum f](k y , r) 
in an implicit way, via r. Such a self-similar dependence on k x is manifested, at each fixed 
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FIG. 3. Spectrum narrowing for large time scales. 

k y , in shrinking of the original k x profile along the /c^-axis as time grows (see Fig. [3]). The 
spectrum is narrowing and its derivative is growing near small values of k x , and when it is 
so steep that a significant variation occurs over the range ~ k*, the kinetic equation breaks 
down. Time estimate for such a breakdown is t ~ {k*)~ 2 - 

Therefore, the weak turbulence description will break down at the late evolution stages, 
and the wave-kinetic equation will no longer work. However, it is possible to amend this 
description to take into account the strongly nonlinear effects and develop a qualitative 
theory of subsequent evolution leading to a steady state. Below we will present a qualitative 
argument which will allow us to obtain such a theory. 

First of all we note that the three- wave interaction is never exactly resonant: it involves 
all the quasi-resonant within a certain small distance from the exact resonant frequency - 
so called nonlinear resonance broadening T ~ t~l . In the other words, the delta function 
in kinetic equation 8{2k\ x ) = 5 (u;& — oj\ — u 2 ) should be replaced by a peaked function 
f(ki x ) with a small but finite width T. For sufficiently smooth spectra the difference from 
the delta-function can be ignored, but for sharp and narrow spectra the integrand in the 
kinetic equation's integral become itself peaked and the delta-function broadening becomes 
important. Its main effect is that of a filter f(k x ) in the k x variable which acts to smoothen 
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any sharp changes over the range Ak ~ k*. The energy is no longer conserved separately at 
each fixed k x . 

For large wavenumbers k x 3> k* (where the spectrum remains slowly varying even 
when it is steep at k x ~ k*) the kinetic equation could be easily amended by replacing 
n T (0,k y ,0) = f n T (k lx , k ly ,r) 5(k lx )dk lx with (n T )(k y ,r) = f n T (k lx ,k ly ,r) f(k lx )dki x . 
For small wavenumbers, k x ~ k*, the effect of the resonance broadening is not reduced to 
such a simple modification of just one function in the integral. It is clear that the spectrum 
at k x = 0, which was fixed in the wave-kinetic approximation, will suffer changes caused by 
the smoothing in the direction determined by spectral slope at small k x : if the gradient is 
positive (negative) the value will increase (decrease), as illustrated in Fig. HI The details of 
the evolution at the small wavenumbers are not important because the combined action of 
the self-similar shrinking and smoothing will lead to a very rapid wipeout of all the gradients 
in k x and formation of a steady state with 77 independent of k x . Correspondingly, the values 
of r) at k x = will adjust themselves to the values at k x = 00. After this moment, when the 
rapid dependencies on k x disappear, the kinetic equation in its usual weak turbulence form 
is valid once again, and it can be used for finding the final steady state spectrum. Since 77 is 
now independent of k x , the steady state could be readily obtained from the formal condition 
i](k y , 0) = t](k y , 00) which simply means our solution independent of r; it has nothing to do 
with the initial/final values of the spectrum in time or k x . 

Thus, the evolution can be summarized as follows. At the early stage, t <C (k*)~ 2 , the 
evolution is described by the three-wave kinetic equation. Then at the advanced stage, 
with characteristic times scales t ~ (fc*) -2 , the kinetic equation is broken down by its own 
evolution. Smoothing of strong gradients in k x occurs, which results in spectrum stabilisation 
and re-emergence of the kinetic equation description at the large time scales t 3> (k*)~ 2 . 
This kinetic equation describes spectrum evolution within the steady state regime. Let us 
now consider the properties of such a steady state. 

A. Spectrum evolution in the steady-state 

Let us now analyse the steady state. Based on what was said in the end of the previous 
section, we will seek a r-independent solution of the pseudo- Fourier space equation (1271) : 

£ 2 (y) - £(y) (8(0) + a d ) + T(y) = (44) 
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FIG. 4. Gradients smoothing process. Four iterations for spectrum value stabilisation are presented on 
this figure. At the first stage, the gradient of spectrum in the vicinity of the k x = (see point 1) is 
positive, the initial value of spectrum r](k y , 0) increases and reaches point 2, and then, after crossing the 
maximum, it moves to point 3, which corresponds to negative slope of the spectrum. Then, the initial 
value decreases and arrives at the position 4. This process will continue until the spectrum stabilises at 
T)(ky,0) = rj(ky, oo). 

(formally coinciding with the condition £(y,0) = £(y,oo) = £{y))- Considering this equa- 
tion at y = we have £(0) = .F(0)/cr<i. Also we have £(0) = J r](ky)dky > (because 
7]{k y ) > 0). Solving the quadratic equation, we have: 



1 — 1 / ^ ^ N 1/2 

£( y ) = ~{a d + ?{oy<T d ) ± - (k + mi^f - my)) 



(45) 



To satisfy condition £ (0) = j^O)/^, we must choose "+" if ^(0) > a\ and "-" otherwise. 

We suppose that the forcing decays at infinity, lim^oo J-{y) — > 0, which is the case eg. 
for the Gaussian forcing. This means that we do not force the k y = mode. Then we see 
that \im y _> 00 £(y) — > if J'(O) < a\ and Xmi y ^ 00 £{y) — > a d + ^(O)/^ if ^(0) > a\ . In 
the second case we have a spectrum with a delta function at k y = 0. Thus we observe an 
interesting phenomenon of condensation into k y = mode in the cases when the forcing 
prevails over the dissipation at the small scales y (corresponding to high k y s). In the first 
case £(y) is a monotonously decreasing (to 0) function of y, where as in the second case it 
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is monotonously increasing (asymptoting to constant). 

The first case is physically more relevant, because in most cases of interest dissipation 
dominates over forcing at the small scales. In this case £(y) behaves qualitatively similar 
as in the two examples considered in section IIVI Namely, if we take the same Gaussian 
forcing as in these two examples, we will have £(y) which has a maximum at y — 0, smooth 
everywhere (including y — 0) and rapidly decaying to zero for y — > oo. However, there is a 
big difference from the previous examples in that now the characteristic width of function 
£(y), and respectively the width of the spectrum in the k y variable, is of the same order as 
the width of the forcing function. Therefore, there is no inertial range in the final steady 
state considered here. This is an even stronger case of the nonlocal interaction than in the 
two examples considered before. Both the forcing and the dissipation parameters enter in 
the final answer, but not the parameters of the initial condition: the steady state beyond 
the weak turbulence has already forgotten all the initial data. 

VI. SUMMARY. 

We have shown that the three-wave interactions for the pseudo-Alfven waves (PAW) in 
the 2D MHD system are non-empty and it is possible to obtain a three-wave kinetic equation 
within the weak turbulence approach. These interactions take place in the second order of 
the anisotropy parameter. 

We found Kolmogorov-Zakharov power law spectra for PAW in 2D MHD system and 
showed that they are not realisable due to divergence of the collision integrals of the kinetic 
equation. In the balanced case this divergence is marginal. This is an indirect indication 
that the 2D PAW turbulence is nonlocal: it is dominated by interaction of waves with very 
different in size wavelengths. Our full analytical solution of the kinetic equation confirms 
such a nonlocality. It also dispels the myth that all marginally nonlocal spectra can be 
"fixed" by a logarithmic correction. 

The crucial technique for our analysis is passing to pseudo-physical space via Fourier- 
transforming the kinetic equation and using a self-similar effective time variable. This has 
allowed us to dramatically simplify the kinetic equation, solve it analytically in some im- 
portant cases and to fully analyse in the other important cases. The two main examples we 
analysed have a Gaussian-shaped forcing of low wavenumbers and a dissipation represented 
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by either uniform friction or by viscosity. The first case is solvable analytically, and the 
second one is shown to possess a similar behaviour. Namely, the spectrum evolves indepen- 
dently at each k x and it tends to a flat steady-state spectrum in the inertial range (which is 
not a log-corrected Kolmogorov-Zakharov spectrum). 

At each fixed k y , the spectrum develops sharp gradients in k x at small k x , which even- 
tually leads to the breakdown of the weak turbulence description. We present a qualitative 
argument about what follows after this moment. We argue that the effect of strong tur- 
bulence is to smoothen the sharp gradients via the nonlinear resonance broadening effect. 
This will lead to a steady state with no gradient in k x for which the weak turbulence kinetic 
equation formally works once again, and we present an analytical solution for such a steady 
state. 

On the practical side, one should derive from our work a warning that the 2D and the 
3D MHD systems are dramatically different, and one should be careful when extrapolating 
the 2D results, eg. numerical ones, onto the 3D case. Indeed, in contrast to 2D, in 3D there 
is no gradient sharpening at small parallel wave numbers, and the Kolmogorov-Zakharov 
spectrum is a local and well behaved solution. 
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Appendix A: Derivation of the wave-kinetic equation 

In section Hi Al we wrote the 2D MHD system in the interaction representation (|T3|) . which 
comprises a starting point 

for derivation of the wave-kinetic equation. Let us define the wave spectrum as 

where the average is taken over the random initial conditions, and L 2 is the area of the 
periodic box. With this normalisation nk tends to a finite limit ~ e 2 as L — > oo provided 
that the wave density is finite and uniform in the 2D physical space. 
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The next step consists in making use of the time scales separation. We are introducing 
an intermediate time scale T which should be much smaller than the typical non-linear 
time, t n i = 27r/(e 2 a;), and much greater than the linear wave period, = 2tt/u>. Taking 
T = 2n/{euj) will satisfy these conditions, ^ < T C t„|. Then we are looking for solutions 
at time t = T in the form series in small e: 



c-(r)=4 ± '°> + 4 ± ' l » + e 2 4 ± - 2 ) + ..., 



(Al) 



where we suppose that the lowest order amplitudes cjjp ^ = c^(0) correspond to the linear 
regime. 

For the spectrum we have: 



[nt(T)-nf(0)]/(e^ 



+ c 



,(±,0)* (±,2) 
J k c k 



(±,0) (±,2)« 
c k c k 



After substituting expansion (lAip into equation (fTH]) in the first order we have: 



.(±,i) 



(T) = XV Mt2 A T (±2A; lie ) cf ' 0) 4 ±,0) 5i 2) 



1,2 



where 



A r (±2fc la 



3 ±i2k lx T _ j 



±2zA; 



1.x 



For the second order we can write: 



(A2) 



(A3) 



(A4) 



V W 8 v V 2t3A 6^' 0) c^ 0) c { t' 0) E(±2k lx ,±2k 3x ) 



1,2,3,4 



+Vi,3,444 ±,0) 4 ±,C M T ' 0) £ (±2k lx , T 2k 3x ) 



(A5) 



with 



E{x,y) 



e* xt A t (y)dt. 



(A6) 



Next, we are going to assume that the initial amplitudes c^' ^ are Gaussian random 
variables which are statistically independent at each k, and use Wick's rule: 



(±,0) (±,0) ( T ,0) (T,0) 
c l c 2 c 3 c 4 



tf(k 1 + ka)*(k B + k4)<|ci ± ' 0) | 2 ><|cS :F ' 0) | a > 



(A7) 



We also remember that because the physical space amplitudes are real functions we have 
(c( ± '°)(k))* = c( ± '°)(-k). 
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We have: 



,(±,i) 



V kth2 V kA4 A T (±2k lx )A* T (±2k 3x ) 



1,2,3,4 



X ( C 



„(=F,0) J±,0) / (T,0) 



-1 



,(±.o) 



°12°34 



(A8) 



12 



and 



(±,0)V f±,2) 



1234 



^,1,2*12 ["^,3,4^34^ (±2*1*, ±2fc 3a; j ^ ^C fc 

O |^4,i, 2 | 2 £ (±2fci*, T2A:i 2 ) nfnf^ 



(±,0)\ =F,0 =F,0 ±,0 



c l c 3 C 4 



(A9) 
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where we have used abbreviations nf = n T (ki, t), = ^ T (k 2 , t) and = 5(ki + k 2 — k). 
Next we note that 



SE(±2k lx ,T2k lx ) = -%E{T2k lx ,±2k 



l.r J 



and 



^(±2^, ^2^) = sRE(=F2fci a: , ±2fc lx ) 



sin 2 (fc^T) 
2 k lx 



Let substitute expressions ( 1A8I) and ( 1A9j) into the eq. ( 1A2I) . 



(T) - nf (0) = T2 I^WI V (4 - O 5 



fc sin (k lx T) 

12" 



1,2 



^la 



(A10) 



(All) 



(A12) 



where we have used that | A T (2fc lx )| 2 = sin 2 (fc la .T)//c 2 3 ,. 

Now we take the infinite box limit, L — > oo, and pass to the continuous description in 
the k-space using the rule 



J2^2 6 ^ ~* I 5 12^k!dk 2 , 

1,2 



where 5\ 2 in the integrand means Dirac's delta (recall that it is Kronecker delta in the sum). 

At the next stage of the wave-kinetic procedure we need to use the weakness of the non- 
linearity in our system by taking the limit e — > 0, which is equivalent to T — > oo. For the 
r.h.s. of the eq. (1A12I) . we obtain: 

sin 2 (k lx T) 



lim 

T^oo k 2 



71 



T5(k lx ) 



(A13) 



l.T 
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Then after multiplying both parts of the eq. ( 1A12j) by 1/T, its l.h.s. becomes: 

- am, (Ai4) 

where we took into account that time T is much less than the nonlinear time at which the 
spectrum evolves. 

After these steps we can finally write down the kinetic equation: 



n u = 7r 



J V knKi H - ni] 6 (k - kx - k 2 ) 5(2 k^dkt dk 2 . (A15) 



Appendix B: Locality study for Kolmogorov-Zakharov solutions. 

In order to explore realisability of Kolmogorov-Zakharov spectra, r](k y \ oo)^ oc k v y , we 
need to proceed with a convergence study of the collisional integrals: 

/oo 
5 (k ly + k 2 y ~ ky) ^y^ {^y^ ~ \ky\^)dk ly dk 2 y (Bl) 
-OO 

Let us consider the first one, choosing a + at the exponent of There are three 

singular points: 

1. kly, k 2 y -> OO, 

2. k\ y y 0, k 2 y y k y , 

3. k 2y -> 0, kiy ->■ k y . 

At the first point, we should use the fact that the integral L \x\ u dx converges when 
v < — 1. After substituting k 2y — k y — ki y , we have : 



/oo 
|^ly| + \\k y k±y\ \ky\ ) dk\y. 



Then two cases are possible: 

• when ct_ > 0, the main contribution is made by: 



\k\y\ ++ dk\y, 



which is convergent for a + + ct_ < —1, 
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when a_ < the expression for the main contribution is made by: 

| ^~\.y I I ^y I dfh~^y . 

ji 

and is convergent for a + < —1. 

At the second singular point, after integration out k 2y using the 5-function in fIBip . we have: 

/ \hy\ a+ (\k y -k ly \ a - -\k y \ a -)dk ly ~ I k^ +1 dk ly , (B2) 
Jo Jo 

and we get the convergence condition a + > —2. To get this condition, we have performed 
the series expansion: \k y — ki y \ a - = \k y \ a ~(l + a_k ly /k y ) + . . . , and we have used the fact 
that the integral j^x u dx is convergent for v > — 1. 

To obtain the convergence condition for the last singular point, we integrate out k\ y using 
the 5-function: 

/ \k y - k 2y \ a + (\k 2y \ a - -\k y \ a -)dk 2y ~ (B3) 
Jo 

I \k y \ + \k 2y \ dk 2y / \k y \ + "^~ dk 2y . 
Jo Jo 

The second integral is always convergent, and the first one is convergent for > —1. 
Finally, the convergence region for the first collisional integral (1B1|) in the space of indices 

is: 

{{(a + + a_ < -1) n (a_ > 0)} U {(«+ < -1) n (a_ < -0)}} n (a+ > -2) n (a_ > -1) . 

It is represented by the grey trapezoid in Fig. [5j 

To find the convergence zone for the second integral of fIBll) (with a_ in the exponent of 
ki y ) we just take reflection of the convergence zone for the first integral with respect to the 
line ct_ = a + . Finally, to get the convergence conditions for both collisional integrals one 
should take the intersection of the both zones. 

As we can see in Fig. 0such an intersection produces a zero set. There are no power law 
exponents a~ and a + for which both collision integrals would be convergent, and there is a 
single point which corresponds to marginal (logarithmic) divergence, a~ = a + = — 1. This 
point corresponds to the Kolmogorov-Zakharov spectrum in the balanced turbulence case. 
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FIG. 5. Locality study for Kolmogorov-Zakharov spectrum. 



Common wisdom 



481 ] is that such marginally nonlocal spectra can be fixed by a logarithmic 



correction. However, in the main text of this paper we show that this is not the case. 
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